In [1]:
%matplotlib inline
from glob import glob
from soma import aims
import pandas as pd
import numpy as np
from scipy import stats
import os.path as osp
from nilearn import plotting
from IPython.display import display_html, display_jpeg, display
Loading the data
In [2]:
datadir = '/home/grg/upf'
controls = sorted(glob(osp.join(datadir, 'Data/Control/s12*.nii')))
diazepam = sorted(glob(osp.join(datadir, 'Data/Diazepam/s12*.nii')))
rois = glob(osp.join(datadir, 'Data/ROI/*.nii'))
roinames = [each.split('/')[-1][:-4] for each in rois]
In [30]:
data = pd.read_excel(osp.join(datadir, 'Dades.xlsx'), header=1)
table_cont = []
table_diaz = []
# For controls first
for i, cont in enumerate(controls):
line = []
c = aims.read(cont).arraydata() # Loads the subject's image
for j, roiname in enumerate(roinames):
roi = aims.read(rois[j]).arraydata() # Loads the ROI mask
line.append(c[roi==1].mean()) # Gets the mean value over the ROI
table_cont.append(line)
# For diazepam subjects then
for i, diaz in enumerate(diazepam):
line = []
d = aims.read(diaz).arraydata()
for j, roiname in enumerate(roinames):
roi = aims.read(rois[j]).arraydata()
line.append(d[roi==1].mean())
table_diaz.append(line)
In [29]:
pd.DataFrame(table_cont, columns=roinames).head()
Out[29]:
In [26]:
table = table_cont[:]
table.extend(table_diaz)
print len(table_cont), len(table_diaz)
roivalues = pd.DataFrame(table, columns=roinames)
groups = ['controls'] * 31
groups.extend(['patients'] * 31)
roivalues['group'] = groups
roivalues.groupby('group').head()
Out[26]:
Converting values to SUV...
In [8]:
def convert(val, p, d):
cf = 65.98
voxel_size = .008
nci_to_bq = 37.0
val = val * cf / nci_to_bq / voxel_size
val = val / (d / p)
return val
In [39]:
ttests = []
suv_controls = []
suv_diaz = []
suv_controls_wb = []
suv_controls_wm = []
# Converting kg to g and mCi to nCi
weights_cont = data['Pes [Kg]'] * 1000.0
dose_cont = data['Dosi [mCi]'] * 1000000.0
weights_diaz = data['Pes [Kg].1'] * 1000.0
dose_diaz = data['Dosi [mCi].1'] * 1000000.0
for i, h in enumerate(roinames):
# Converting values to SUV for controls
a = np.array(np.array(table_cont)[:,i], dtype=np.float64)
a = convert(a, weights_cont, dose_cont)
suv_controls.append(a)
# Converting values to SUV for diazepam
b = np.array(np.array(table_diaz)[:,i], dtype=np.float64)
b = convert(b, weights_diaz, dose_diaz)
suv_diaz.append(b)
# Computing t-tests
ttests.append(stats.ttest_ind(a, b)[1] / 2.0)
In [40]:
suv_controls = np.array(suv_controls).transpose()
suv_diaz = np.array(suv_diaz).transpose()
Now what effects do we observe ?
In [43]:
pd.DataFrame(ttests, index=roinames, columns=['suv'])
Out[43]:
In [44]:
def plot_data(suv_controls, suv_diaz, ylim=None):
from matplotlib import pyplot as plt
import matplotlib
matplotlib.rc('figure', facecolor='white')
fig, ax1 = plt.subplots(figsize=(10, 6))
pos_cont = xrange(0,len(roinames)*2,2)
pos_diaz = xrange(1,len(roinames)*2,2)
box_c = plt.boxplot(np.array(suv_controls), positions=pos_cont, patch_artist=True)
for patch in box_c['boxes']:
patch.set(facecolor='cyan')
box_d = plt.boxplot(np.array(suv_diaz), positions=pos_diaz, patch_artist=True)
for patch in box_d['boxes']:
patch.set(facecolor='pink')
ax1.yaxis.grid(True, linestyle='-', which='major', color='lightgray',
alpha=0.5)
ax1.set_xlim(-1)
if not ylim is None:
ax1.set_ylim([0,ylim])
plt.xticks(rotation=70)
ax1.set_xticks(list(np.array(range(0, 2*len(roinames), 2)) + np.array([0.5]*len(roinames))))
ax1.set_xticklabels(roinames)
plt.figtext(0.80, 0.0015, 'controls', backgroundcolor='cyan', color='black', weight='roman', size='x-small')
plt.figtext(0.80, 0.025, 'diazepam', backgroundcolor='pink', color='black', weight='roman', size='x-small')
plt.show()
In [45]:
plot_data(suv_controls, suv_diaz)
In [49]:
from matplotlib import pyplot as plt
_ = plt.boxplot(np.array(suv_controls))
In [ ]: